function F = cal_chg_ge(eqm, firm, eqm_cf, firm_cf)
%% Counterfactual changes by ex-ante quintile

    % Indicators for ex-ante quintile
    Qindicator = eqm.Qindicator; %Nf by 5 matrix
    
    
    % Useful functions
    my_fcn1 = @(x, density) sum(kron(ones(1,5),x.*density).*Qindicator,1)./sum(kron(ones(1,5),density).*Qindicator,1)*100;
    my_fcn2 = @(x, density) sum(kron(ones(1,5),x.*density).*kron(ones(2,1),Qindicator),1)./sum(kron(ones(1,5),density).*kron(ones(2,1),Qindicator),1)*100; 
    my_fcn3 = @(x, density) sum(kron(ones(1,5),x.*density).*kron(ones(3,1),Qindicator),1)./sum(kron(ones(1,5),density).*kron(ones(3,1),Qindicator),1)*100; 

    % Quality
    logdiff = log(firm_cf.qlevel) - log(firm.qlevel);  
    ge_cf.eqt_qlevel_switcher      = my_fcn1(logdiff, eqm_cf.density_switcher);
    ge_cf.eqt_qlevel_exp_ctn       = my_fcn1(logdiff, eqm_cf.density_exp_ctn);
    ge_cf.eqt_qlevel_exp           = my_fcn1(logdiff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);
    ge_cf.eqt_qlevel_nonexp_ctn    = my_fcn1(logdiff, eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_qlevel_all           = my_fcn1(logdiff, eqm_cf.density_all);

    % Productivity
    logdiff = firm_cf.lnz - firm.lnz;  
    ge_cf.eqt_z_switcher      = my_fcn1(logdiff, eqm_cf.density_switcher);
    ge_cf.eqt_z_exp_ctn       = my_fcn1(logdiff, eqm_cf.density_exp_ctn);
    ge_cf.eqt_z_exp           = my_fcn1(logdiff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);
    ge_cf.eqt_z_nonexp_ctn    = my_fcn1(logdiff, eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_z_all           = my_fcn1(logdiff, eqm_cf.density_all);   
    
    % Wage
    logdiff = log(firm_cf.wage) - log(firm.wage);   
    ge_cf.eqt_wage_switcher    = my_fcn1(logdiff, eqm_cf.density_switcher);
    ge_cf.eqt_wage_exp_ctn     = my_fcn1(logdiff, eqm_cf.density_exp_ctn);
    ge_cf.eqt_wage_exp         = my_fcn1(logdiff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);  
    ge_cf.eqt_wage_nonexp_ctn  = my_fcn1(logdiff, eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_wage_all         = my_fcn1(logdiff, eqm_cf.density_all);
    
    % Sales
    logdiff_exp_ctn     = firm_cf.lnsales_exp - firm.lnsales_exp;
    logdiff_nonexp_ctn  = firm_cf.lnsales_nonexp - firm.lnsales_nonexp;
    logdiff_switcher    = firm_cf.lnsales_exp - firm.lnsales_nonexp;  
    ge_cf.eqt_sales_switcher   = my_fcn1(logdiff_switcher,    eqm_cf.density_switcher);
    ge_cf.eqt_sales_exp_ctn    = my_fcn1(logdiff_exp_ctn,     eqm_cf.density_exp_ctn);
    ge_cf.eqt_sales_exp        = my_fcn2([logdiff_exp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_switcher]);
    ge_cf.eqt_sales_nonexp_ctn = my_fcn1(logdiff_nonexp_ctn,  eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_sales_all        = my_fcn3([logdiff_exp_ctn; logdiff_nonexp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_nonexp_ctn; eqm_cf.density_switcher]);
    
    % Number of suppliers
    logdiff_exp_ctn     = log(firm_cf.indegree_exp) - log(firm.indegree_exp);
    logdiff_nonexp_ctn  = log(firm_cf.indegree_nonexp) - log(firm.indegree_nonexp);
    logdiff_switcher    = log(firm_cf.indegree_exp) - log(firm.indegree_nonexp);   
    ge_cf.eqt_indegree_switcher   = my_fcn1(logdiff_switcher,     eqm_cf.density_switcher);
    ge_cf.eqt_indegree_exp_ctn    = my_fcn1(logdiff_exp_ctn,      eqm_cf.density_exp_ctn);
    ge_cf.eqt_indegree_exp        = my_fcn2([logdiff_exp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_switcher]);
    ge_cf.eqt_indegree_nonexp_ctn = my_fcn1(logdiff_nonexp_ctn,   eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_indegree_all        = my_fcn3([logdiff_exp_ctn; logdiff_nonexp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_nonexp_ctn; eqm_cf.density_switcher]);

    % Number of customers
    logdiff_exp_ctn     = log(firm_cf.outdegree_exp) - log(firm.outdegree_exp);
    logdiff_nonexp_ctn  = log(firm_cf.outdegree_nonexp) - log(firm.outdegree_nonexp);
    logdiff_switcher    = log(firm_cf.outdegree_exp) - log(firm.outdegree_nonexp); 
    ge_cf.eqt_outdegree_switcher   = my_fcn1(logdiff_switcher,    eqm_cf.density_switcher);
    ge_cf.eqt_outdegree_exp_ctn    = my_fcn1(logdiff_exp_ctn,     eqm_cf.density_exp_ctn);
    ge_cf.eqt_outdegree_exp        = my_fcn2([logdiff_exp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_switcher]);
    ge_cf.eqt_outdegree_nonexp_ctn = my_fcn1(logdiff_nonexp_ctn,  eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_outdegree_all        = my_fcn3([logdiff_exp_ctn; logdiff_nonexp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_nonexp_ctn; eqm_cf.density_switcher]);
        
    % Number of workers
    logdiff_exp_ctn     = log(firm_cf.employment_exp) - log(firm.employment_exp);
    logdiff_nonexp_ctn  = log(firm_cf.employment_nonexp) - log(firm.employment_nonexp);
    logdiff_switcher    = log(firm_cf.employment_exp) - log(firm.employment_nonexp);  
    ge_cf.eqt_employment_switcher   = my_fcn1(logdiff_switcher,   eqm_cf.density_switcher);
    ge_cf.eqt_employment_exp_ctn    = my_fcn1(logdiff_exp_ctn,    eqm_cf.density_exp_ctn);
    ge_cf.eqt_employment_exp        = my_fcn2([logdiff_exp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_switcher]);
    ge_cf.eqt_employment_nonexp_ctn = my_fcn1(logdiff_nonexp_ctn, eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_employment_all        = my_fcn3([logdiff_exp_ctn; logdiff_nonexp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_nonexp_ctn; eqm_cf.density_switcher]);
    
       
    % Average log wage of suppliers
    logdiff = firm_cf.avg_wgt_lnwageS - firm.avg_wgt_lnwageS;
    ge_cf.eqt_lnwageS_switcher      = my_fcn1(logdiff, eqm_cf.density_switcher);
    ge_cf.eqt_lnwageS_exp_ctn       = my_fcn1(logdiff, eqm_cf.density_exp_ctn);
    ge_cf.eqt_lnwageS_exp           = my_fcn1(logdiff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);
    ge_cf.eqt_lnwageS_nonexp_ctn    = my_fcn1(logdiff, eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_lnwageS_all           = my_fcn1(logdiff, eqm_cf.density_all);    
    
    
    % Fraction of higher quality suppliers
    logdiff = log(firm_cf.frac_S_higher_q) - log(firm.frac_S_higher_q);   
    ge_cf.eqt_frac_S_higher_switcher      = my_fcn1(logdiff, eqm_cf.density_switcher);
    ge_cf.eqt_frac_S_higher_exp_ctn       = my_fcn1(logdiff, eqm_cf.density_exp_ctn);
    ge_cf.eqt_frac_S_higher_exp           = my_fcn1(logdiff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);
    ge_cf.eqt_frac_S_higher_nonexp_ctn    = my_fcn1(logdiff, eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_frac_S_higher_all           = my_fcn1(logdiff, eqm_cf.density_all);    
    
    
    % Fraction of higher quality customers
    logdiff = log(firm_cf.frac_B_higher_q) - log(firm.frac_B_higher_q);   
    ge_cf.eqt_frac_B_higher_switcher      = my_fcn1(logdiff, eqm_cf.density_switcher);
    ge_cf.eqt_frac_B_higher_exp_ctn       = my_fcn1(logdiff, eqm_cf.density_exp_ctn);
    ge_cf.eqt_frac_B_higher_exp           = my_fcn1(logdiff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);
    ge_cf.eqt_frac_B_higher_nonexp_ctn    = my_fcn1(logdiff, eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_frac_B_higher_all           = my_fcn1(logdiff, eqm_cf.density_all);    
    


    % Export probability
    diff = firm_cf.prob_exp - firm.prob_exp;  
    ge_cf.eqt_exprob_all           = my_fcn1(diff, eqm_cf.density_all);      
    
        % July27 check
        ge_cf.check_exp_prob_cf = my_fcn1(firm_cf.prob_exp, eqm_cf.density_all);
        ge_cf.check_exp_prob_init = my_fcn1(firm.prob_exp, eqm_cf.density_all);

    
    % Export intensity
    diff = firm_cf.exp_intst_exp - firm.exp_intst_exp;  
    ge_cf.eqt_exintst_exp          = my_fcn1(diff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);
   
        % July27 check
        ge_cf.check_exp_intst_init = my_fcn1(firm.exp_intst_exp, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);
        ge_cf.check_exp_intst_cf = my_fcn1(firm_cf.exp_intst_exp, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);


    
    
%% Overall counterfactual changes
       
    % Useful function
    my_avg = @(x, density) sum(x.*density)./sum(density);

    % Quality
    ratio = firm_cf.qlevel./firm.qlevel;
    ge_cf.overall_qlevel_switcher   = (my_avg(ratio, eqm_cf.density_switcher)-1)*100;
    ge_cf.overall_qlevel_exp_ctn    = (my_avg(ratio, eqm_cf.density_exp_ctn)-1)*100;
    ge_cf.overall_qlevel_exp        = (my_avg(ratio, eqm_cf.density_exp_ctn+eqm_cf.density_switcher)-1)*100;
    ge_cf.overall_qlevel_nonexp_ctn = (my_avg(ratio, eqm_cf.density_nonexp_ctn)-1)*100;
    ge_cf.overall_qlevel_all        = (my_avg(ratio, eqm_cf.density_all)-1)*100;
       
    % Wage
    ratio = firm_cf.wage./firm.wage;
    ge_cf.overall_wage_switcher   = (my_avg(ratio, eqm_cf.density_switcher)-1)*100;
    ge_cf.overall_wage_exp_ctn    = (my_avg(ratio, eqm_cf.density_exp_ctn)-1)*100;
    ge_cf.overall_wage_exp        = (my_avg(ratio, eqm_cf.density_exp_ctn+eqm_cf.density_switcher)-1)*100;
    ge_cf.overall_wage_nonexp_ctn = (my_avg(ratio, eqm_cf.density_nonexp_ctn)-1)*100;
    ge_cf.overall_wage_all        = (my_avg(ratio, eqm_cf.density_all)-1)*100;    
    
    % Average log wage of suppliers
    logdiff = firm_cf.avg_wgt_lnwageS - firm.avg_wgt_lnwageS;
    ge_cf.overall_lnwageS_switcher      = my_avg(logdiff, eqm_cf.density_switcher)*100;
    ge_cf.overall_lnwageS_exp_ctn       = my_avg(logdiff, eqm_cf.density_exp_ctn)*100;
    ge_cf.overall_lnwageS_exp           = my_avg(logdiff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher)*100;
    ge_cf.overall_lnwageS_nonexp_ctn    = my_avg(logdiff, eqm_cf.density_nonexp_ctn)*100;
    ge_cf.overall_lnwageS_all           = my_avg(logdiff, eqm_cf.density_all)*100;
    



%% Break down changes in output (%)

    ge_cf.sales_total   = (eqm_cf.X - eqm.X)*100;
    ge_cf.sales_manuf   = (sum(eqm_cf.X_m)-sum(eqm.X_m))*100;
    ge_cf.sales_foreign = (eqm_cf.B1-eqm.B1)*100;
    ge_cf.sales_service = ( sum(eqm_cf.D_s./eqm_cf.D0.*eqm_cf.X0+eqm_cf.D_s.*eqm_cf.rv_ads./eqm_cf.D1.*eqm_cf.X1) ...
                           -sum(eqm.D_s./eqm.D0.*eqm.X0+eqm.D_s.*eqm.rv_ads./eqm.D1.*eqm.X1) )*100;



%% Return

    F = ge_cf;

    
end